##Replication code for "Understanding Political Communication and Political Communicators on Twitch"
##Reference network and description data object creation
##Sangyeon Kim

library(igraph)

#1.top20 network creation

load("./raw_chatpost_data/top20_bychannel.rda")

#dropping out channels with only one chat, that are all not in the list
top20_bychannel[[22]]<-NULL
top20_bychannel[[21]]<-NULL
top20_bychannel[[20]]<-NULL
top20_bychannel[[19]]<-NULL

networkinfo<-vector(mode="list",length=18)

for(i in 1:18){
  networkinfo[[i]]<-vector(mode="list",length=5)
}

#create edge_chats
for(i in 1:length(networkinfo)){
  a<-top20_bychannel[[i]]
  channelname<-a$channelname
  networkinfo[[i]][[1]]<-unique(channelname)
  b<-grepl('@',a$text)
  
  c<-which(b==TRUE)
  
  nodelist<-unique(a$username)
  
  d<-a[c,]
  
  #exclude those from bots
  
  e<-grepl('bot',nodelist)
  
  f<-which(e==TRUE)
  if(length(f)>0){
    g<-nodelist[f]
    
    h<-d[!d$username==g[1],]
    for(k in 2:length(g)){
      h<-h[!h$username==g[k],]
    }
    edge_chats<-h
    rm(a,b,c,d,e,f,g,h)
  }
  else{
    edge_chats<-d
    rm(a,b,c,d,e,f)
  }
  
  
  
  
  ######
  
  #create sender, receiver, and text objects
  
  library(stringr)
  sender<-vector(mode="character",length=length(edge_chats$chatid))
  receiver<-vector(mode="list",length=length(edge_chats$chatid))
  text<-vector(mode="character",length=length(edge_chats$chatid))
  for(p in 1:length(sender)){
    a<-edge_chats[p,]
    sender[[p]]<-a$username
    b<-str_split(a$text," ")
    c<-str_detect(b[[1]],"@")
    d<-b[[1]][which(c==TRUE)]
    e<-b[[1]][which(c==FALSE)]
    receiver[[p]]<-tolower(str_replace_all(d, "@", ""))
    text[[p]]<-str_flatten(e, collapse = " ")
    print(p)
  }
  rm(a,b,c,d,e,p)
  
  #create edgelist
  edgelist_uag<-vector(mode="list",length=length(sender))
  for(q in 1:length(edgelist_uag)){
    
    a<-sender[q]
    b<-receiver[q]
    c<-text[q]
    d<-cbind(a,b[[1]][1],c)
    if(length(b[[1]])>2){
      for(t in 2:length(b[[1]])){
        d<-rbind(d,cbind(a,b[[1]][t],c))
      }
    }
    edgelist_uag[[q]]<-d
    print(q)
  }
  
  a<-edgelist_uag[[1]]
  for(u in 2:length(edgelist_uag)){
    a<-rbind(a,edgelist_uag[[u]])
  }
  
  edgelist_wtext<-a
  rm(a,b,c,k)
  colnames(edgelist_wtext)<-c("sender","receiver","text")
  
  #sort out edgelist with receiver that are not in the nodelist
  edgelist_wtext<-as.data.frame(edgelist_wtext)
  
  edgelist_wtext$sender<-as.character(edgelist_wtext$sender)
  edgelist_wtext$receiver<-as.character(edgelist_wtext$receiver)
  
  h<-edgelist_wtext[,1] %in% nodelist
  h2<-edgelist_wtext[,2] %in% nodelist 
  length(which(h)) #28776
  length(which(h2)) #27948
  
  edgelist_wtext<-edgelist_wtext[which(h2),]
  h<-edgelist_wtext[,1] %in% nodelist
  h2<-edgelist_wtext[,2] %in% nodelist 
  length(which(h)) #27948
  length(which(h2)) #27948 #now it all matches!
  
  #create edgelist without text
  edgelist<-edgelist_wtext[,c(1,2)]
  
  
  networkinfo[[i]][[2]]<-nodelist
  networkinfo[[i]][[3]]<-edgelist
  networkinfo[[i]][[4]]<-edgelist_wtext
  
  #creating network
  
  ####layout function to be used for plotting####
  
  library(igraph)
  g <- graph.data.frame(edgelist, nodelist, directed=T)
  networkinfo[[i]][[5]]<-g
  print(i)
}

save(networkinfo,file="./network_data/top20_networkinfo.rda")

#2.top21 to top100 network creation

load("./raw_chatpost_data/top20_100_bychannel.rda")

networkinfo_2<-vector(mode="list",length=length(top20_100_bychannel))

for(i in 1:length(networkinfo_2)){
  networkinfo_2[[i]]<-vector(mode="list",length=5)
}

#create edge_chats
for(i in 1:length(networkinfo_2)){
  a<-top20_100_bychannel[[i]]
  channelname<-a$channelname
  networkinfo_2[[i]][[1]]<-channelname
  b<-grepl('@',a$text)
  
  c<-which(b==TRUE)
  
  nodelist<-unique(a$username)
  
  d<-a[c,]
  
  #exclude those from bots
  
  e<-grepl('bot',nodelist)
  
  f<-which(e==TRUE)
  if(length(f)>0){
    g<-nodelist[f]
    
    h<-d[!d$username==g[1],]
    for(k in 2:length(g)){
      h<-h[!h$username==g[k],]
    }
    edge_chats<-h
    rm(a,b,c,d,e,f,g,h)
  }
  else{
    edge_chats<-d
    rm(a,b,c,d,e,f)
  }
  
  
  
  
  ######
  
  #create sender, receiver, and text objects
  
  library(stringr)
  sender<-vector(mode="character",length=length(edge_chats$chatid))
  receiver<-vector(mode="list",length=length(edge_chats$chatid))
  text<-vector(mode="character",length=length(edge_chats$chatid))
  for(p in 1:length(sender)){
    a<-edge_chats[p,]
    sender[[p]]<-a$username
    b<-str_split(a$text," ")
    c<-str_detect(b[[1]],"@")
    d<-b[[1]][which(c==TRUE)]
    e<-b[[1]][which(c==FALSE)]
    receiver[[p]]<-tolower(str_replace_all(d, "@", ""))
    text[[p]]<-str_flatten(e, collapse = " ")
  }
  rm(a,b,c,d,e,p)
  
  #create edgelist
  edgelist_uag<-vector(mode="list",length=length(sender))
  for(q in 1:length(edgelist_uag)){
    
    a<-sender[q]
    b<-receiver[q]
    c<-text[q]
    d<-cbind(a,b[[1]][1],c)
    if(length(b[[1]])>2){
      for(t in 2:length(b[[1]])){
        d<-rbind(d,cbind(a,b[[1]][t],c))
      }
    }
    edgelist_uag[[q]]<-d
    
  }
  
  a<-edgelist_uag[[1]]
  for(u in 2:length(edgelist_uag)){
    a<-rbind(a,edgelist_uag[[u]])
  }
  
  edgelist_wtext<-a
  rm(a,b,c,k)
  colnames(edgelist_wtext)<-c("sender","receiver","text")
  
  #sort out edgelist with receiver that are not in the nodelist
  edgelist_wtext<-as.data.frame(edgelist_wtext)
  
  edgelist_wtext$sender<-as.character(edgelist_wtext$sender)
  edgelist_wtext$receiver<-as.character(edgelist_wtext$receiver)
  
  h<-edgelist_wtext[,1] %in% nodelist
  h2<-edgelist_wtext[,2] %in% nodelist 
  length(which(h)) #28776
  length(which(h2)) #27948
  
  edgelist_wtext<-edgelist_wtext[which(h2),]
  h<-edgelist_wtext[,1] %in% nodelist
  h2<-edgelist_wtext[,2] %in% nodelist 
  length(which(h)) #27948
  length(which(h2)) #27948 #now it all matches!
  
  #create edgelist without text
  edgelist<-edgelist_wtext[,c(1,2)]
  
  
  networkinfo_2[[i]][[2]]<-nodelist
  networkinfo_2[[i]][[3]]<-edgelist
  networkinfo_2[[i]][[4]]<-edgelist_wtext
  
  #creating network
  
  ####layout function to be used for plotting####
  
  library(igraph)
  g <- graph.data.frame(edgelist, nodelist, directed=T)
  networkinfo_2[[i]][[5]]<-g
  print(i)
}

save(networkinfo_2,file="./network_data/top21_100_networkinfo.rda")


#3. top101 to top337 network creation

load("./raw_chatpost_data/top101_337_bychannel.rda")

#before go in, get rid of the data with only one observation, as it does not show any meaningful patterns
oneornot<-vector(mode="logical",length=length(top101_337_bychannel))
for(i in 1:length(oneornot)){
  oneornot[[i]]<-length(top101_337_bychannel[[i]]$username)>1
}
one<-which(oneornot==FALSE)

top101_337_bychannel[[one[[4]]]]<-NULL
top101_337_bychannel[[one[[3]]]]<-NULL
top101_337_bychannel[[one[[2]]]]<-NULL
top101_337_bychannel[[one[[1]]]]<-NULL

networkinfo_3<-vector(mode="list",length=length(top101_337_bychannel))

for(i in 1:length(networkinfo_3)){
  networkinfo_3[[i]]<-vector(mode="list",length=5)
}


#create edge_chats
for(i in 1:length(networkinfo_3)){
  a<-top101_337_bychannel[[i]]
  channelname<-a$channelname
  networkinfo_3[[i]][[1]]<-unique(channelname)
  b<-grepl('@',a$text)
  
  c<-which(b==TRUE)
  
  nodelist<-unique(a$username)
  
  d<-a[c,]
  
  #exclude those from bots
  
  e<-grepl('bot',nodelist)
  
  f<-which(e==TRUE)
  if(length(f)>0){
    g<-nodelist[f]
    
    h<-d[!d$username==g[1],]
    for(k in 2:length(g)){
      h<-h[!h$username==g[k],]
    }
    edge_chats<-h
    rm(a,b,c,d,e,f,g,h)
  }
  else{
    edge_chats<-d
    rm(a,b,c,d,e,f)
  }
  
  
  
  
  ######
  
  #create sender, receiver, and text objects
  
  library(stringr)
  sender<-vector(mode="character",length=length(edge_chats$chatid))
  receiver<-vector(mode="list",length=length(edge_chats$chatid))
  text<-vector(mode="character",length=length(edge_chats$chatid))
  if(length(text)>0){
    for(p in 1:length(sender)){
      a<-edge_chats[p,]
      sender[[p]]<-a$username
      b<-str_split(a$text," ")
      c<-str_detect(b[[1]],"@")
      d<-b[[1]][which(c==TRUE)]
      e<-b[[1]][which(c==FALSE)]
      receiver[[p]]<-tolower(str_replace_all(d, "@", ""))
      text[[p]]<-str_flatten(e, collapse = " ")
    }
    rm(a,b,c,d,e,p)
  }
  
  #create edgelist
  if(length(text)>0){
    edgelist_uag<-vector(mode="list",length=length(sender))
    for(q in 1:length(edgelist_uag)){
      
      a<-sender[q]
      b<-receiver[q]
      c<-text[q]
      d<-cbind(a,b[[1]][1],c)
      if(length(b[[1]])>2){
        for(t in 2:length(b[[1]])){
          d<-rbind(d,cbind(a,b[[1]][t],c))
        }
      }
      edgelist_uag[[q]]<-d
      
    }
    
    a<-edgelist_uag[[1]]
    if(length(edgelist_uag)>1){
      for(u in 2:length(edgelist_uag)){
        a<-rbind(a,edgelist_uag[[u]])
      }
    }
    edgelist_wtext<-a
    rm(a,b,c,k)
    colnames(edgelist_wtext)<-c("sender","receiver","text")
    
    #sort out edgelist with receiver that are not in the nodelist
    edgelist_wtext<-as.data.frame(edgelist_wtext)
    
    edgelist_wtext$sender<-as.character(edgelist_wtext$sender)
    edgelist_wtext$receiver<-as.character(edgelist_wtext$receiver)
    
    h<-edgelist_wtext[,1] %in% nodelist
    h2<-edgelist_wtext[,2] %in% nodelist 
    length(which(h)) #28776
    length(which(h2)) #27948
    
    edgelist_wtext<-edgelist_wtext[which(h2),]
    h<-edgelist_wtext[,1] %in% nodelist
    h2<-edgelist_wtext[,2] %in% nodelist 
    length(which(h)) #27948
    length(which(h2)) #27948 #now it all matches!
    
    #create edgelist without text
    edgelist<-edgelist_wtext[,c(1,2)]
    
    
    networkinfo_3[[i]][[2]]<-nodelist
    networkinfo_3[[i]][[3]]<-edgelist
    networkinfo_3[[i]][[4]]<-edgelist_wtext
    
    #creating network
    
    ####layout function to be used for plotting####
    
    library(igraph)
    g <- graph.data.frame(edgelist, nodelist, directed=T)
    networkinfo_3[[i]][[5]]<-g
  }
  else{
    edgelist <- data.frame(sender=character(),
                           receiver=character())
    networkinfo_3[[i]][[2]]<-nodelist
    networkinfo_3[[i]][[3]]<-edgelist
    g <- graph.data.frame(edgelist, nodelist, directed=T)
    networkinfo_3[[i]][[5]]<-g
  }
  print(i)
}
save(networkinfo_3,file="./network_data/top101_337_networkinfo.rda")

#4.top338 to top574

load("./raw_chatpost_data/top338_574_bychannel.rda")
#before go in, get rid of the data with only one observation, as it does not show any meaningful patterns
oneornot<-vector(mode="logical",length=length(top338_574_bychannel))
for(i in 1:length(oneornot)){
  oneornot[[i]]<-length(top338_574_bychannel[[i]]$username)>1
}
one<-which(oneornot==FALSE)
top338_574_bychannel[[one[[8]]]]<-NULL
top338_574_bychannel[[one[[7]]]]<-NULL
top338_574_bychannel[[one[[6]]]]<-NULL
top338_574_bychannel[[one[[5]]]]<-NULL
top338_574_bychannel[[one[[4]]]]<-NULL
top338_574_bychannel[[one[[3]]]]<-NULL
top338_574_bychannel[[one[[2]]]]<-NULL
top338_574_bychannel[[one[[1]]]]<-NULL

networkinfo_4<-vector(mode="list",length=length(top338_574_bychannel))

for(i in 1:length(networkinfo_4)){
  networkinfo_4[[i]]<-vector(mode="list",length=5)
}

#create edge_chats
for(i in 1:length(networkinfo_4)){
  a<-top338_574_bychannel[[i]]
  channelname<-a$channelname
  networkinfo_4[[i]][[1]]<-unique(channelname)
  b<-grepl('@',a$text)
  
  c<-which(b==TRUE)
  
  nodelist<-unique(a$username)
  
  d<-a[c,]
  
  #exclude those from bots
  
  e<-grepl('bot',nodelist)
  
  f<-which(e==TRUE)
  if(length(f)>0){
    g<-nodelist[f]
    
    h<-d[!d$username==g[1],]
    for(k in 2:length(g)){
      h<-h[!h$username==g[k],]
    }
    edge_chats<-h
    rm(a,b,c,d,e,f,g,h)
  }
  else{
    edge_chats<-d
    rm(a,b,c,d,e,f)
  }
  
  
  
  
  ######
  
  #create sender, receiver, and text objects
  
  library(stringr)
  sender<-vector(mode="character",length=length(edge_chats$chatid))
  receiver<-vector(mode="list",length=length(edge_chats$chatid))
  text<-vector(mode="character",length=length(edge_chats$chatid))
  if(length(text)>0){
    for(p in 1:length(sender)){
      a<-edge_chats[p,]
      sender[[p]]<-a$username
      b<-str_split(a$text," ")
      c<-str_detect(b[[1]],"@")
      d<-b[[1]][which(c==TRUE)]
      e<-b[[1]][which(c==FALSE)]
      receiver[[p]]<-tolower(str_replace_all(d, "@", ""))
      text[[p]]<-str_flatten(e, collapse = " ")
    }
    rm(a,b,c,d,e,p)
  }
  
  #create edgelist
  if(length(text)>0){
    edgelist_uag<-vector(mode="list",length=length(sender))
    for(q in 1:length(edgelist_uag)){
      
      a<-sender[q]
      b<-receiver[q]
      c<-text[q]
      d<-cbind(a,b[[1]][1],c)
      if(length(b[[1]])>2){
        for(t in 2:length(b[[1]])){
          d<-rbind(d,cbind(a,b[[1]][t],c))
        }
      }
      edgelist_uag[[q]]<-d
      
    }
    
    a<-edgelist_uag[[1]]
    if(length(edgelist_uag)>1){
      for(u in 2:length(edgelist_uag)){
        a<-rbind(a,edgelist_uag[[u]])
      }
    }
    edgelist_wtext<-a
    rm(a,b,c,k)
    colnames(edgelist_wtext)<-c("sender","receiver","text")
    
    #sort out edgelist with receiver that are not in the nodelist
    edgelist_wtext<-as.data.frame(edgelist_wtext)
    
    edgelist_wtext$sender<-as.character(edgelist_wtext$sender)
    edgelist_wtext$receiver<-as.character(edgelist_wtext$receiver)
    
    h<-edgelist_wtext[,1] %in% nodelist
    h2<-edgelist_wtext[,2] %in% nodelist 
    length(which(h)) #28776
    length(which(h2)) #27948
    
    edgelist_wtext<-edgelist_wtext[which(h2),]
    h<-edgelist_wtext[,1] %in% nodelist
    h2<-edgelist_wtext[,2] %in% nodelist 
    length(which(h)) #27948
    length(which(h2)) #27948 #now it all matches!
    
    #create edgelist without text
    edgelist<-edgelist_wtext[,c(1,2)]
    
    
    networkinfo_4[[i]][[2]]<-nodelist
    networkinfo_4[[i]][[3]]<-edgelist
    networkinfo_4[[i]][[4]]<-edgelist_wtext
    
    #creating network
    
    ####layout function to be used for plotting####
    
    library(igraph)
    g <- graph.data.frame(edgelist, nodelist, directed=T)
    networkinfo_4[[i]][[5]]<-g
  }
  else{
    edgelist <- data.frame(sender=character(),
                           receiver=character())
    networkinfo_4[[i]][[2]]<-nodelist
    networkinfo_4[[i]][[3]]<-edgelist
    g <- graph.data.frame(edgelist, nodelist, directed=T)
    networkinfo_4[[i]][[5]]<-g
  }
  print(i)
}

save(networkinfo_4,file="./network_data/top338_574_networkinfo.rda")

#5. Network description stat object creation

load("./network_data/top20_networkinfo.rda")
load("./network_data/top21_100_networkinfo.rda")
load("./network_data/top101_337_networkinfo.rda")
load("./network_data/top338_574_networkinfo.rda")

network_sumstat<-vector(mode="list",length=length(networkinfo)+length(networkinfo_2)+length(networkinfo_3)+length(networkinfo_4))

for(i in 1:length(network_sumstat)){
  network_sumstat[[i]]<-vector(mode="list",length=14)
}

list3<-c(3,4,5,7,8,9,10,11,12)

for(i in 1:length(network_sumstat)){
  for(j in 1:length(list3)){
    network_sumstat[[i]][[list3[j]]]<-vector(mode="list",length=3)
  }
}

for(i in 2:length(networkinfo)){
  print(i)
  g<-networkinfo[[i]][[5]]
  
  #centrality: degree centrality
  channel<-unique(networkinfo[[i]][[1]])
  network_sumstat[[i]][[1]]<-unique(networkinfo[[i]][[1]])
  node<-networkinfo[[i]][[2]]
  
  owner_chat<-length(which(channel==node)) #see whether the channel owner participate in streaming chat 1=yes, 0=no
  network_sumstat[[i]][[2]]<-owner_chat
  
  if(length(networkinfo[[i]][[3]]$sender)>0){
    deg_all<-igraph::degree(g)
    network_sumstat[[i]][[3]][[1]]<-deg_all
    deg_all_his<-hist(deg_all, xlab="Degree",main="Degree Distribution (All)", col="Gray")
    network_sumstat[[i]][[3]][[2]]<-deg_all_his
    deg_all_log_his<-hist(log(deg_all), xlab="Degree (logged)",main="Degree Distribution (All - logged)", col="Gray")
    network_sumstat[[i]][[3]][[3]]<-deg_all_log_his
    
    deg_out<-igraph::degree(g, mode="out")
    network_sumstat[[i]][[4]][[1]]<-deg_out
    deg_out_his<-hist(deg_out, xlab="Degree (out)",main="Degree Distribution (Out)", col="Gray")
    network_sumstat[[i]][[4]][[2]]<-deg_out_his
    deg_out_log_his<-hist(log(deg_out), xlab="Degree (logged)",main="Degree Distribution (All - logged)", col="Gray")
    network_sumstat[[i]][[4]][[3]]<-deg_out_log_his
    
    deg_in<-igraph::degree(g, mode="in")
    network_sumstat[[i]][[5]][[1]]<-deg_in
    deg_in_his<-hist(deg_in, xlab="Degree (in)",main="Degree Distribution (In)", col="Gray")
    network_sumstat[[i]][[5]][[2]]<-deg_in_his
    deg_in_log_his<-hist(log(deg_in), xlab="Degree (logged)",main="Degree Distribution (In - logged)", col="Gray")
    network_sumstat[[i]][[5]][[3]]<-deg_in_log_his
    #reciprocity
    
    recip_ratio<-reciprocity(g,mode='ratio')
    network_sumstat[[i]][[6]]<-recip_ratio
    
    #popularity
    
    #1) Modularity/community memebership
    # measures the strength of division of a network into clusters and cummunities (modules).
    # Networks with high modularity have dense connections between the nodes within modules 
    # but sparse connections between nodes in different modules. 
    # The “fast greedy” method, starts with each node belonging to a separate community 
    # and these are itteratively merged such that it yields the largest increase in 
    # the current value of modularity.
    
    # Find communities using walktrap algorithm 
    community_walktrap <- walktrap.community(g, 
                                             steps = 4, modularity=TRUE,membership=TRUE) # The length of the random walks to perform.
    network_sumstat[[i]][[7]][[1]]<-community_walktrap
    # list of nodes in the most densely connected subgraph extracted from the community_walktrap object
    membership <- membership(community_walktrap)
    network_sumstat[[i]][[7]][[2]]<-membership
    modularity<- modularity(community_walktrap)
    network_sumstat[[i]][[7]][[3]]<-modularity
    #Modularity scores of +1 mean that all the edges in a community are connecting nodes within the community.
    #A score of 0 would mean that the community has half its edges connecting nodes within the same community,
    #and half connecting nodes outside the community.
    #A score of -1 means that there are no edges connecting nodes within the community,
    # and they instead all connect nodes outside the community.
    
    ################################
    #graph that shows community structure of the network
    #for(i in 1:length(a[[1]])){
    #  if(a[[1]][[i]]>5){
    #    a[[1]][[i]]<-NA
    #  }
    #  else{}
    #}
    
    #b<-a[[1]]
    #b<-na.omit(b)
    #a[[1]]<-b
    #mycol <- t_col("blue", perc = 98.5, name = "lt.blue")
    #V(g)$label <- NA
    #V(g)$size<-0.1
    #
    #E(g)$color<-mycol 
    #a<-vcount(g)^2
    #large<-layout_with_lgl(
    #  g,
    #  maxiter = 150,
    #  maxdelta = 5000,
    #  area = a,
    #  coolexp = 1.5,
    #  repulserad = a * 5000,
    #  cellsize = sqrt(sqrt(a)),
    #  root = NULL
    #)
    
    #jpeg("test.jpg", width = 1200, height = 1200)
    #plot(g, 
    #     mark.groups=a, #A list of numeric vectors. The communities can be highlighted using colored polygons.      
    #     mark.col="pink",
    #     layout=large)
    #dev.off()
    #################################
    
    #2) log-log plot to see whether the network follows power law
    
    llplot_all<-plot(1:length(deg_all_his$counts),deg_all_his$counts+1, xlab = "Degree (all)", ylab = "Frequency", cex.lab = 1.5, main = "Degree-all (log-log scale)", log = "xy", type = "b")
    network_sumstat[[i]][[8]][[1]]<-llplot_all
    llplot_out<-plot(1:length(deg_out_his$counts),deg_out_his$counts+1, xlab = "Degree (out)", ylab = "Frequency", cex.lab = 1.5, main = "Degree-out (log-log scale)", log = "xy", type = "b")
    network_sumstat[[i]][[8]][[2]]<-llplot_out
    llplot_in<-plot(1:length(deg_in_his$counts),deg_in_his$counts+1, xlab = "Degree (out)", ylab = "Frequency", cex.lab = 1.5, main = "Degree-in (log-log scale)", log = "xy", type = "b")
    network_sumstat[[i]][[8]][[3]]<-llplot_in
    #function to test whether the graph shows power law or not
    #?fit_power_law
    power_all <- fit_power_law(deg_all)
    network_sumstat[[i]][[9]][[1]]<-power_all
    power_out <- fit_power_law(deg_out)
    network_sumstat[[i]][[9]][[2]]<-power_out
    power_in <- fit_power_law(deg_in)
    network_sumstat[[i]][[9]][[3]]<-power_in
    
    # exponent/alpha
    power_all$alpha
    network_sumstat[[i]][[10]][[1]]<-power_all$alpha
    power_out$alpha
    network_sumstat[[i]][[10]][[2]]<-power_out$alpha
    power_in$alpha
    network_sumstat[[i]][[10]][[3]]<-power_in$alpha
    # test statistic of a Kolmogorov-Smirnov test 
    # Kolmogorov-Smirnov test. Remember if p value is less than 0.05, 
    #we can reject the null that our network data could have been drawn from a fitted power-law distribution.
    network_sumstat[[i]][[11]][[1]]<-power_all$KS.stat
    network_sumstat[[i]][[11]][[2]]<-power_out$KS.stat
    network_sumstat[[i]][[11]][[3]]<-power_in$KS.stat
    # p value for the K.S. test at 5%
    network_sumstat[[i]][[12]][[1]]<-power_all$KS.p
    network_sumstat[[i]][[12]][[2]]<-power_out$KS.p
    network_sumstat[[i]][[12]][[3]]<-power_in$KS.p
    
    #transitivity
    trans<-transitivity(g,isolates = "zero") 
    network_sumstat[[i]][[13]]<-trans
    #$T = 1$ implies perfect transitivity, i.e., a network whose components are all cliques.
    #$T=0$ implied no closed path of length two, which happens for various topologies, such as a tree or a square lattice.
    #assortative mixing
    deg_hom<-assortativity_degree(g)
    #over 0, means there is homophily among nodes with same/similar degree
    network_sumstat[[i]][[14]]<-deg_hom
  }
}


for(i in 1:length(networkinfo_2)){
  print(i)
  g<-networkinfo_2[[i]][[5]]
  
  #centrality: degree centrality
  channel<-unique(networkinfo_2[[i]][[1]])
  network_sumstat[[i+length(networkinfo)]][[1]]<-unique(networkinfo_2[[i]][[1]])
  node<-networkinfo_2[[i]][[2]]
  
  owner_chat<-length(which(channel==node)) #see whether the channel owner participate in streaming chat 1=yes, 0=no
  network_sumstat[[i+length(networkinfo)]][[2]]<-owner_chat
  
  if(length(networkinfo_2[[i]][[3]]$sender)>0){
    deg_all<-igraph::degree(g)
    network_sumstat[[i+length(networkinfo)]][[3]][[1]]<-deg_all
    deg_all_his<-hist(deg_all, xlab="Degree",main="Degree Distribution (All)", col="Gray")
    network_sumstat[[i+length(networkinfo)]][[3]][[2]]<-deg_all_his
    deg_all_log_his<-hist(log(deg_all), xlab="Degree (logged)",main="Degree Distribution (All - logged)", col="Gray")
    network_sumstat[[i+length(networkinfo)]][[3]][[3]]<-deg_all_log_his
    
    deg_out<-igraph::degree(g, mode="out")
    network_sumstat[[i+length(networkinfo)]][[4]][[1]]<-deg_out
    deg_out_his<-hist(deg_out, xlab="Degree (out)",main="Degree Distribution (Out)", col="Gray")
    network_sumstat[[i+length(networkinfo)]][[4]][[2]]<-deg_out_his
    deg_out_log_his<-hist(log(deg_out), xlab="Degree (logged)",main="Degree Distribution (All - logged)", col="Gray")
    network_sumstat[[i+length(networkinfo)]][[4]][[3]]<-deg_out_log_his
    
    deg_in<-igraph::degree(g, mode="in")
    network_sumstat[[i+length(networkinfo)]][[5]][[1]]<-deg_in
    deg_in_his<-hist(deg_in, xlab="Degree (in)",main="Degree Distribution (In)", col="Gray")
    network_sumstat[[i+length(networkinfo)]][[5]][[2]]<-deg_in_his
    deg_in_log_his<-hist(log(deg_in), xlab="Degree (logged)",main="Degree Distribution (In - logged)", col="Gray")
    network_sumstat[[i+length(networkinfo)]][[5]][[3]]<-deg_in_log_his
    #reciprocity
    
    recip_ratio<-reciprocity(g,mode='ratio')
    network_sumstat[[i+length(networkinfo)]][[6]]<-recip_ratio
    
    #popularity
    
    #1) Modularity/community memebership
    # measures the strength of division of a network into clusters and cummunities (modules).
    # Networks with high modularity have dense connections between the nodes within modules 
    # but sparse connections between nodes in different modules. 
    # The “fast greedy” method, starts with each node belonging to a separate community 
    # and these are itteratively merged such that it yields the largest increase in 
    # the current value of modularity.
    
    # Find communities using walktrap algorithm 
    community_walktrap <- walktrap.community(g, 
                                             steps = 4, modularity=TRUE,membership=TRUE) # The length of the random walks to perform.
    network_sumstat[[i+length(networkinfo)]][[7]][[1]]<-community_walktrap
    # list of nodes in the most densely connected subgraph extracted from the community_walktrap object
    membership <- membership(community_walktrap)
    network_sumstat[[i+length(networkinfo)]][[7]][[2]]<-membership
    modularity<- modularity(community_walktrap)
    network_sumstat[[i+length(networkinfo)]][[7]][[3]]<-modularity
    #Modularity scores of +1 mean that all the edges in a community are connecting nodes within the community.
    #A score of 0 would mean that the community has half its edges connecting nodes within the same community,
    #and half connecting nodes outside the community.
    #A score of -1 means that there are no edges connecting nodes within the community,
    # and they instead all connect nodes outside the community.
    
    ################################
    #graph that shows community structure of the network
    #for(i in 1:length(a[[1]])){
    #  if(a[[1]][[i]]>5){
    #    a[[1]][[i]]<-NA
    #  }
    #  else{}
    #}
    
    #b<-a[[1]]
    #b<-na.omit(b)
    #a[[1]]<-b
    #mycol <- t_col("blue", perc = 98.5, name = "lt.blue")
    #V(g)$label <- NA
    #V(g)$size<-0.1
    #
    #E(g)$color<-mycol 
    #a<-vcount(g)^2
    #large<-layout_with_lgl(
    #  g,
    #  maxiter = 150,
    #  maxdelta = 5000,
    #  area = a,
    #  coolexp = 1.5,
    #  repulserad = a * 5000,
    #  cellsize = sqrt(sqrt(a)),
    #  root = NULL
    #)
    
    #jpeg("test.jpg", width = 1200, height = 1200)
    #plot(g, 
    #     mark.groups=a, #A list of numeric vectors. The communities can be highlighted using colored polygons.      
    #     mark.col="pink",
    #     layout=large)
    #dev.off()
    #################################
    
    #2) log-log plot to see whether the network follows power law
    
    llplot_all<-plot(1:length(deg_all_his$counts),deg_all_his$counts+1, xlab = "Degree (all)", ylab = "Frequency", cex.lab = 1.5, main = "Degree-all (log-log scale)", log = "xy", type = "b")
    network_sumstat[[i+length(networkinfo)]][[8]][[1]]<-llplot_all
    llplot_out<-plot(1:length(deg_out_his$counts),deg_out_his$counts+1, xlab = "Degree (out)", ylab = "Frequency", cex.lab = 1.5, main = "Degree-out (log-log scale)", log = "xy", type = "b")
    network_sumstat[[i+length(networkinfo)]][[8]][[2]]<-llplot_out
    llplot_in<-plot(1:length(deg_in_his$counts),deg_in_his$counts+1, xlab = "Degree (out)", ylab = "Frequency", cex.lab = 1.5, main = "Degree-in (log-log scale)", log = "xy", type = "b")
    network_sumstat[[i+length(networkinfo)]][[8]][[3]]<-llplot_in
    #function to test whether the graph shows power law or not
    #?fit_power_law
    power_all <- fit_power_law(deg_all)
    network_sumstat[[i+length(networkinfo)]][[9]][[1]]<-power_all
    power_out <- fit_power_law(deg_out)
    network_sumstat[[i+length(networkinfo)]][[9]][[2]]<-power_out
    power_in <- fit_power_law(deg_in)
    network_sumstat[[i+length(networkinfo)]][[9]][[3]]<-power_in
    
    # exponent/alpha
    power_all$alpha
    network_sumstat[[i+length(networkinfo)]][[10]][[1]]<-power_all$alpha
    power_out$alpha
    network_sumstat[[i+length(networkinfo)]][[10]][[2]]<-power_out$alpha
    power_in$alpha
    network_sumstat[[i+length(networkinfo)]][[10]][[3]]<-power_in$alpha
    # test statistic of a Kolmogorov-Smirnov test 
    # Kolmogorov-Smirnov test. Remember if p value is less than 0.05, 
    #we can reject the null that our network data could have been drawn from a fitted power-law distribution.
    network_sumstat[[i+length(networkinfo)]][[11]][[1]]<-power_all$KS.stat
    network_sumstat[[i+length(networkinfo)]][[11]][[2]]<-power_out$KS.stat
    network_sumstat[[i+length(networkinfo)]][[11]][[3]]<-power_in$KS.stat
    # p value for the K.S. test at 5%
    network_sumstat[[i+length(networkinfo)]][[12]][[1]]<-power_all$KS.p
    network_sumstat[[i+length(networkinfo)]][[12]][[2]]<-power_out$KS.p
    network_sumstat[[i+length(networkinfo)]][[12]][[3]]<-power_in$KS.p
    
    #transitivity
    trans<-transitivity(g,isolates = "zero") 
    network_sumstat[[i+length(networkinfo)]][[13]]<-trans
    #$T = 1$ implies perfect transitivity, i.e., a network whose components are all cliques.
    #$T=0$ implied no closed path of length two, which happens for various topologies, such as a tree or a square lattice.
    #assortative mixing
    deg_hom<-assortativity_degree(g)
    #over 0, means there is homophily among nodes with same/similar degree
    network_sumstat[[i+length(networkinfo)]][[14]]<-deg_hom
  }
}

for(i in 1:length(networkinfo_3)){
  print(i)
  g<-networkinfo_3[[i]][[5]]
  
  #centrality: degree centrality
  channel<-unique(networkinfo_3[[i]][[1]])
  network_sumstat[[i+length(networkinfo)+length(networkinfo_2)]][[1]]<-unique(networkinfo_3[[i]][[1]])
  node<-networkinfo_3[[i]][[2]]
  
  owner_chat<-length(which(channel==node)) #see whether the channel owner participate in streaming chat 1=yes, 0=no
  network_sumstat[[i+length(networkinfo)+length(networkinfo_2)]][[2]]<-owner_chat
  
  if(length(networkinfo_3[[i]][[3]]$sender)>0){
    deg_all<-igraph::degree(g)
    network_sumstat[[i+length(networkinfo)+length(networkinfo_2)]][[3]][[1]]<-deg_all
    deg_all_his<-hist(deg_all, xlab="Degree",main="Degree Distribution (All)", col="Gray")
    network_sumstat[[i+length(networkinfo)+length(networkinfo_2)]][[3]][[2]]<-deg_all_his
    deg_all_log_his<-hist(log(deg_all), xlab="Degree (logged)",main="Degree Distribution (All - logged)", col="Gray")
    network_sumstat[[i+length(networkinfo)+length(networkinfo_2)]][[3]][[3]]<-deg_all_log_his
    
    deg_out<-igraph::degree(g, mode="out")
    network_sumstat[[i+length(networkinfo)+length(networkinfo_2)]][[4]][[1]]<-deg_out
    deg_out_his<-hist(deg_out, xlab="Degree (out)",main="Degree Distribution (Out)", col="Gray")
    network_sumstat[[i+length(networkinfo)+length(networkinfo_2)]][[4]][[2]]<-deg_out_his
    deg_out_log_his<-hist(log(deg_out), xlab="Degree (logged)",main="Degree Distribution (All - logged)", col="Gray")
    network_sumstat[[i+length(networkinfo)+length(networkinfo_2)]][[4]][[3]]<-deg_out_log_his
    
    deg_in<-igraph::degree(g, mode="in")
    network_sumstat[[i+length(networkinfo)+length(networkinfo_2)]][[5]][[1]]<-deg_in
    deg_in_his<-hist(deg_in, xlab="Degree (in)",main="Degree Distribution (In)", col="Gray")
    network_sumstat[[i+length(networkinfo)+length(networkinfo_2)]][[5]][[2]]<-deg_in_his
    deg_in_log_his<-hist(log(deg_in), xlab="Degree (logged)",main="Degree Distribution (In - logged)", col="Gray")
    network_sumstat[[i+length(networkinfo)+length(networkinfo_2)]][[5]][[3]]<-deg_in_log_his
    #reciprocity
    
    recip_ratio<-reciprocity(g,mode='ratio')
    network_sumstat[[i+length(networkinfo)+length(networkinfo_2)]][[6]]<-recip_ratio
    
    #popularity
    
    #1) Modularity/community memebership
    # measures the strength of division of a network into clusters and cummunities (modules).
    # Networks with high modularity have dense connections between the nodes within modules 
    # but sparse connections between nodes in different modules. 
    # The “fast greedy” method, starts with each node belonging to a separate community 
    # and these are itteratively merged such that it yields the largest increase in 
    # the current value of modularity.
    
    # Find communities using walktrap algorithm 
    community_walktrap <- walktrap.community(g, 
                                             steps = 4, modularity=TRUE,membership=TRUE) # The length of the random walks to perform.
    network_sumstat[[i+length(networkinfo)+length(networkinfo_2)]][[7]][[1]]<-community_walktrap
    # list of nodes in the most densely connected subgraph extracted from the community_walktrap object
    membership <- membership(community_walktrap)
    network_sumstat[[i+length(networkinfo)+length(networkinfo_2)]][[7]][[2]]<-membership
    modularity<- modularity(community_walktrap)
    network_sumstat[[i+length(networkinfo)+length(networkinfo_2)]][[7]][[3]]<-modularity
    #Modularity scores of +1 mean that all the edges in a community are connecting nodes within the community.
    #A score of 0 would mean that the community has half its edges connecting nodes within the same community,
    #and half connecting nodes outside the community.
    #A score of -1 means that there are no edges connecting nodes within the community,
    # and they instead all connect nodes outside the community.
    
    ################################
    #graph that shows community structure of the network
    #for(i in 1:length(a[[1]])){
    #  if(a[[1]][[i]]>5){
    #    a[[1]][[i]]<-NA
    #  }
    #  else{}
    #}
    
    #b<-a[[1]]
    #b<-na.omit(b)
    #a[[1]]<-b
    #mycol <- t_col("blue", perc = 98.5, name = "lt.blue")
    #V(g)$label <- NA
    #V(g)$size<-0.1
    #
    #E(g)$color<-mycol 
    #a<-vcount(g)^2
    #large<-layout_with_lgl(
    #  g,
    #  maxiter = 150,
    #  maxdelta = 5000,
    #  area = a,
    #  coolexp = 1.5,
    #  repulserad = a * 5000,
    #  cellsize = sqrt(sqrt(a)),
    #  root = NULL
    #)
    
    #jpeg("test.jpg", width = 1200, height = 1200)
    #plot(g, 
    #     mark.groups=a, #A list of numeric vectors. The communities can be highlighted using colored polygons.      
    #     mark.col="pink",
    #     layout=large)
    #dev.off()
    #################################
    
    #2) log-log plot to see whether the network follows power law
    
    llplot_all<-plot(1:length(deg_all_his$counts),deg_all_his$counts+1, xlab = "Degree (all)", ylab = "Frequency", cex.lab = 1.5, main = "Degree-all (log-log scale)", log = "xy", type = "b")
    network_sumstat[[i+length(networkinfo)+length(networkinfo_2)]][[8]][[1]]<-llplot_all
    llplot_out<-plot(1:length(deg_out_his$counts),deg_out_his$counts+1, xlab = "Degree (out)", ylab = "Frequency", cex.lab = 1.5, main = "Degree-out (log-log scale)", log = "xy", type = "b")
    network_sumstat[[i+length(networkinfo)+length(networkinfo_2)]][[8]][[2]]<-llplot_out
    llplot_in<-plot(1:length(deg_in_his$counts),deg_in_his$counts+1, xlab = "Degree (out)", ylab = "Frequency", cex.lab = 1.5, main = "Degree-in (log-log scale)", log = "xy", type = "b")
    network_sumstat[[i+length(networkinfo)+length(networkinfo_2)]][[8]][[3]]<-llplot_in
    #function to test whether the graph shows power law or not
    #?fit_power_law
    power_all <- fit_power_law(deg_all)
    network_sumstat[[i+length(networkinfo)+length(networkinfo_2)]][[9]][[1]]<-power_all
    power_out <- fit_power_law(deg_out)
    network_sumstat[[i+length(networkinfo)+length(networkinfo_2)]][[9]][[2]]<-power_out
    power_in <- fit_power_law(deg_in)
    network_sumstat[[i+length(networkinfo)+length(networkinfo_2)]][[9]][[3]]<-power_in
    
    # exponent/alpha
    power_all$alpha
    network_sumstat[[i+length(networkinfo)+length(networkinfo_2)]][[10]][[1]]<-power_all$alpha
    power_out$alpha
    network_sumstat[[i+length(networkinfo)+length(networkinfo_2)]][[10]][[2]]<-power_out$alpha
    power_in$alpha
    network_sumstat[[i+length(networkinfo)+length(networkinfo_2)]][[10]][[3]]<-power_in$alpha
    # test statistic of a Kolmogorov-Smirnov test 
    # Kolmogorov-Smirnov test. Remember if p value is less than 0.05, 
    #we can reject the null that our network data could have been drawn from a fitted power-law distribution.
    network_sumstat[[i+length(networkinfo)+length(networkinfo_2)]][[11]][[1]]<-power_all$KS.stat
    network_sumstat[[i+length(networkinfo)+length(networkinfo_2)]][[11]][[2]]<-power_out$KS.stat
    network_sumstat[[i+length(networkinfo)+length(networkinfo_2)]][[11]][[3]]<-power_in$KS.stat
    # p value for the K.S. test at 5%
    network_sumstat[[i+length(networkinfo)+length(networkinfo_2)]][[12]][[1]]<-power_all$KS.p
    network_sumstat[[i+length(networkinfo)+length(networkinfo_2)]][[12]][[2]]<-power_out$KS.p
    network_sumstat[[i+length(networkinfo)+length(networkinfo_2)]][[12]][[3]]<-power_in$KS.p
    
    #transitivity
    trans<-transitivity(g,isolates = "zero") 
    network_sumstat[[i+length(networkinfo)+length(networkinfo_2)]][[13]]<-trans
    #$T = 1$ implies perfect transitivity, i.e., a network whose components are all cliques.
    #$T=0$ implied no closed path of length two, which happens for various topologies, such as a tree or a square lattice.
    #assortative mixing
    deg_hom<-assortativity_degree(g)
    #over 0, means there is homophily among nodes with same/similar degree
    network_sumstat[[i+length(networkinfo)+length(networkinfo_2)]][[14]]<-deg_hom
  }
}

for(i in 1:length(networkinfo_4)){
  print(i)
  g<-networkinfo_4[[i]][[5]]
  
  #centrality: degree centrality
  channel<-unique(networkinfo_4[[i]][[1]])
  network_sumstat[[i+length(networkinfo)+length(networkinfo_2)+length(networkinfo_3)]][[1]]<-unique(networkinfo_4[[i]][[1]])
  node<-networkinfo_4[[i]][[2]]
  
  owner_chat<-length(which(channel==node)) #see whether the channel owner participate in streaming chat 1=yes, 0=no
  network_sumstat[[i+length(networkinfo)+length(networkinfo_2)+length(networkinfo_3)]][[2]]<-owner_chat
  
  if(length(networkinfo_4[[i]][[3]]$sender)>0){
    deg_all<-igraph::degree(g)
    network_sumstat[[i+length(networkinfo)+length(networkinfo_2)+length(networkinfo_3)]][[3]][[1]]<-deg_all
    deg_all_his<-hist(deg_all, xlab="Degree",main="Degree Distribution (All)", col="Gray")
    network_sumstat[[i+length(networkinfo)+length(networkinfo_2)+length(networkinfo_3)]][[3]][[2]]<-deg_all_his
    deg_all_log_his<-hist(log(deg_all), xlab="Degree (logged)",main="Degree Distribution (All - logged)", col="Gray")
    network_sumstat[[i+length(networkinfo)+length(networkinfo_2)+length(networkinfo_3)]][[3]][[3]]<-deg_all_log_his
    
    deg_out<-igraph::degree(g, mode="out")
    network_sumstat[[i+length(networkinfo)+length(networkinfo_2)+length(networkinfo_3)]][[4]][[1]]<-deg_out
    deg_out_his<-hist(deg_out, xlab="Degree (out)",main="Degree Distribution (Out)", col="Gray")
    network_sumstat[[i+length(networkinfo)+length(networkinfo_2)+length(networkinfo_3)]][[4]][[2]]<-deg_out_his
    deg_out_log_his<-hist(log(deg_out), xlab="Degree (logged)",main="Degree Distribution (All - logged)", col="Gray")
    network_sumstat[[i+length(networkinfo)+length(networkinfo_2)+length(networkinfo_3)]][[4]][[3]]<-deg_out_log_his
    
    deg_in<-igraph::degree(g, mode="in")
    network_sumstat[[i+length(networkinfo)+length(networkinfo_2)+length(networkinfo_3)]][[5]][[1]]<-deg_in
    deg_in_his<-hist(deg_in, xlab="Degree (in)",main="Degree Distribution (In)", col="Gray")
    network_sumstat[[i+length(networkinfo)+length(networkinfo_2)+length(networkinfo_3)]][[5]][[2]]<-deg_in_his
    deg_in_log_his<-hist(log(deg_in), xlab="Degree (logged)",main="Degree Distribution (In - logged)", col="Gray")
    network_sumstat[[i+length(networkinfo)+length(networkinfo_2)+length(networkinfo_3)]][[5]][[3]]<-deg_in_log_his
    #reciprocity
    
    recip_ratio<-reciprocity(g,mode='ratio')
    network_sumstat[[i+length(networkinfo)+length(networkinfo_2)+length(networkinfo_3)]][[6]]<-recip_ratio
    
    #popularity
    
    #1) Modularity/community memebership
    # measures the strength of division of a network into clusters and cummunities (modules).
    # Networks with high modularity have dense connections between the nodes within modules 
    # but sparse connections between nodes in different modules. 
    # The “fast greedy” method, starts with each node belonging to a separate community 
    # and these are itteratively merged such that it yields the largest increase in 
    # the current value of modularity.
    
    # Find communities using walktrap algorithm 
    community_walktrap <- walktrap.community(g, 
                                             steps = 4, modularity=TRUE,membership=TRUE) # The length of the random walks to perform.
    network_sumstat[[i+length(networkinfo)+length(networkinfo_2)+length(networkinfo_3)]][[7]][[1]]<-community_walktrap
    # list of nodes in the most densely connected subgraph extracted from the community_walktrap object
    membership <- membership(community_walktrap)
    network_sumstat[[i+length(networkinfo)+length(networkinfo_2)+length(networkinfo_3)]][[7]][[2]]<-membership
    modularity<- modularity(community_walktrap)
    network_sumstat[[i+length(networkinfo)+length(networkinfo_2)+length(networkinfo_3)]][[7]][[3]]<-modularity
    #Modularity scores of +1 mean that all the edges in a community are connecting nodes within the community.
    #A score of 0 would mean that the community has half its edges connecting nodes within the same community,
    #and half connecting nodes outside the community.
    #A score of -1 means that there are no edges connecting nodes within the community,
    # and they instead all connect nodes outside the community.
    
    ################################
    #graph that shows community structure of the network
    #for(i in 1:length(a[[1]])){
    #  if(a[[1]][[i]]>5){
    #    a[[1]][[i]]<-NA
    #  }
    #  else{}
    #}
    
    #b<-a[[1]]
    #b<-na.omit(b)
    #a[[1]]<-b
    #mycol <- t_col("blue", perc = 98.5, name = "lt.blue")
    #V(g)$label <- NA
    #V(g)$size<-0.1
    #
    #E(g)$color<-mycol 
    #a<-vcount(g)^2
    #large<-layout_with_lgl(
    #  g,
    #  maxiter = 150,
    #  maxdelta = 5000,
    #  area = a,
    #  coolexp = 1.5,
    #  repulserad = a * 5000,
    #  cellsize = sqrt(sqrt(a)),
    #  root = NULL
    #)
    
    #jpeg("test.jpg", width = 1200, height = 1200)
    #plot(g, 
    #     mark.groups=a, #A list of numeric vectors. The communities can be highlighted using colored polygons.      
    #     mark.col="pink",
    #     layout=large)
    #dev.off()
    #################################
    
    #2) log-log plot to see whether the network follows power law
    
    llplot_all<-plot(1:length(deg_all_his$counts),deg_all_his$counts+1, xlab = "Degree (all)", ylab = "Frequency", cex.lab = 1.5, main = "Degree-all (log-log scale)", log = "xy", type = "b")
    network_sumstat[[i+length(networkinfo)+length(networkinfo_2)+length(networkinfo_3)]][[8]][[1]]<-llplot_all
    llplot_out<-plot(1:length(deg_out_his$counts),deg_out_his$counts+1, xlab = "Degree (out)", ylab = "Frequency", cex.lab = 1.5, main = "Degree-out (log-log scale)", log = "xy", type = "b")
    network_sumstat[[i+length(networkinfo)+length(networkinfo_2)+length(networkinfo_3)]][[8]][[2]]<-llplot_out
    llplot_in<-plot(1:length(deg_in_his$counts),deg_in_his$counts+1, xlab = "Degree (out)", ylab = "Frequency", cex.lab = 1.5, main = "Degree-in (log-log scale)", log = "xy", type = "b")
    network_sumstat[[i+length(networkinfo)+length(networkinfo_2)+length(networkinfo_3)]][[8]][[3]]<-llplot_in
    #function to test whether the graph shows power law or not
    #?fit_power_law
    power_all <- fit_power_law(deg_all)
    network_sumstat[[i+length(networkinfo)+length(networkinfo_2)+length(networkinfo_3)]][[9]][[1]]<-power_all
    power_out <- fit_power_law(deg_out)
    network_sumstat[[i+length(networkinfo)+length(networkinfo_2)+length(networkinfo_3)]][[9]][[2]]<-power_out
    power_in <- fit_power_law(deg_in)
    network_sumstat[[i+length(networkinfo)+length(networkinfo_2)+length(networkinfo_3)]][[9]][[3]]<-power_in
    
    # exponent/alpha
    power_all$alpha
    network_sumstat[[i+length(networkinfo)+length(networkinfo_2)+length(networkinfo_3)]][[10]][[1]]<-power_all$alpha
    power_out$alpha
    network_sumstat[[i+length(networkinfo)+length(networkinfo_2)+length(networkinfo_3)]][[10]][[2]]<-power_out$alpha
    power_in$alpha
    network_sumstat[[i+length(networkinfo)+length(networkinfo_2)+length(networkinfo_3)]][[10]][[3]]<-power_in$alpha
    # test statistic of a Kolmogorov-Smirnov test 
    # Kolmogorov-Smirnov test. Remember if p value is less than 0.05, 
    #we can reject the null that our network data could have been drawn from a fitted power-law distribution.
    network_sumstat[[i+length(networkinfo)+length(networkinfo_2)+length(networkinfo_3)]][[11]][[1]]<-power_all$KS.stat
    network_sumstat[[i+length(networkinfo)+length(networkinfo_2)+length(networkinfo_3)]][[11]][[2]]<-power_out$KS.stat
    network_sumstat[[i+length(networkinfo)+length(networkinfo_2)+length(networkinfo_3)]][[11]][[3]]<-power_in$KS.stat
    # p value for the K.S. test at 5%
    network_sumstat[[i+length(networkinfo)+length(networkinfo_2)+length(networkinfo_3)]][[12]][[1]]<-power_all$KS.p
    network_sumstat[[i+length(networkinfo)+length(networkinfo_2)+length(networkinfo_3)]][[12]][[2]]<-power_out$KS.p
    network_sumstat[[i+length(networkinfo)+length(networkinfo_2)+length(networkinfo_3)]][[12]][[3]]<-power_in$KS.p
    
    #transitivity
    trans<-transitivity(g,isolates = "zero") 
    network_sumstat[[i+length(networkinfo)+length(networkinfo_2)+length(networkinfo_3)]][[13]]<-trans
    #$T = 1$ implies perfect transitivity, i.e., a network whose components are all cliques.
    #$T=0$ implied no closed path of length two, which happens for various topologies, such as a tree or a square lattice.
    #assortative mixing
    deg_hom<-assortativity_degree(g)
    #over 0, means there is homophily among nodes with same/similar degree
    network_sumstat[[i+length(networkinfo)+length(networkinfo_2)+length(networkinfo_3)]][[14]]<-deg_hom
  }
}

save(network_sumstat,file="figure14_to_18.rda")
